/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/operation.hxx"

#ifdef LINAL_DISPATCHER_ENABLE
#include "linal/device/operation.hxx"
#include "linal/device/device.hxx"
#include "linal/device/cache.hxx"
#include "linal/dispatcher.hxx"
#endif

#include "linal/device/host/operation.hxx"

namespace linal {
  /*!
    LU with partial pivoting.
  */
  int lu_piv(Hier_ A, Hier_ p) {
    if (!A.get_m() || !A.get_n()) return true;

    /*
    LINAL_ERROR_DETAIL( A.get_m() == A.get_n(),
                        "%s\n %s %d\n %s %d\n",
                        LINAL_ERROR_A_SQUARE_MATRIX,
                        "length = ", A.get_m(),
                        "width  = ", A.get_n() );
    */

    FLA_Obj ATL,   ATR,      A00,  A01,  A02,
            ABL,   ABR,      A10,  A11,  A12,
                             A20,  A21,  A22;
    FLA_Obj pt, pb,          p0, p1, p2;

    FLA_Part_2x2( ~A, &ATL, &ATR,
                      &ABL, &ABR, 0, 0, FLA_TL );

    FLA_Part_2x1( ~p, &pt,
                      &pb, 0, FLA_TOP );

    while ( FLA_Obj_width( ATL ) < FLA_Obj_width( ~A ) ){
      FLA_Repart_2x2_to_3x3( ATL, ATR,     &A00, &A01, &A02,
                                           &A10, &A11, &A12,
                             ABL, ABR,     &A20, &A21, &A22,
                             1, 1, FLA_BR );
      FLA_Repart_2x1_to_3x1( pt, &p0,
                                 &p1,
                             pb, &p2, 1, FLA_BOTTOM );
      //------------------------------------------------------------
      Hier_
        A_10(A10), A_11(A11), A_12(A12),
        A_21(A21), A_22(A22), p_1(p1);

      FLA_Obj AB0, AB1, AB2;
      FLA_Merge_2x1( A10,
                     A20, &AB0 );
      FLA_Merge_2x1( A11,
                     A21, &AB1 );
      FLA_Merge_2x1( A12,
                     A22, &AB2 );
      Hier_ AB_0(AB0), AB_1(AB1), AB_2(AB2);

#ifdef LINAL_DISPATCHER_ENABLE
      get_dispatcher()->lu_piv( ~(AB_1.flat()), 
                                ((FLA_Obj*)p_1.get_buffer())[0] );
#else
      host::internal::lu_piv( ~(AB_1.flat()), 
                              ((FLA_Obj*)p_1.get_buffer())[0] );
#endif
        
#pragma omp task firstprivate( p_1, AB_0 )
      apply_pivots( FLA_LEFT, FLA_NO_TRANSPOSE, p_1, AB_0 );

#pragma omp task firstprivate( p_1, AB_2 )
      apply_pivots( FLA_LEFT, FLA_NO_TRANSPOSE, p_1, AB_2 );

#pragma omp taskwait

      LINAL_LU_TRSM_A12( A_11, A_12 );
      LINAL_LU_GEMM( A_21, A_12, A_22 );

      //------------------------------------------------------------
      FLA_Cont_with_3x3_to_2x2( &ATL,  &ATR,  A00, A01, A02,
                                              A10, A11, A12,
                                &ABL, &ABR,   A20, A21, A22,
                                FLA_TL );
      FLA_Cont_with_3x1_to_2x1( &pt, p0,
                                     p1,
                                &pb, p2, FLA_TOP );
    }
      
    return true;
  }
}
